Systematic errors due to linear congruential random-number generators 
with the Swendsen— Wang algorithm: A warning 
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We show that hnear congruential pseudo-random-number generators can cause systematic errors 
in Monte Carlo simulations using the Swendsen- Wang algorithm, if the lattice size is a multiple of 
a very large power of 2 and one random number is used per bond. These systematic errors arise 
from correlations within a single bond-update half-sweep. The errors can be eliminated (or at least 
radically reduced) by updating the bonds in a random order or in an aperiodic manner. It also helps 
to use a generator of large modulus (e.g. 60 or more bits) . 
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I. INTRODUCTION 

It has been known for about two decades that linear 
congruential pseudo-random-number generators [ll suffer 
from strong long-range correlations 0, 0, 0, 01^: for 
instance, generators with modulus m = 2^ have strong 
correlations at lags that are multiples of 2*^ whenever 
the ratio k/fi \s large enough Q, y, Furthermore, 
these long-range correlations are known to give rise to 
systematic errors in Monte Carlo simulations employing 
local (e.g. Metropolis or heat bath) updates whenever the 
lattice sites are updated in a fixed order and the number 
of random numbers used per sweep is a multiple of a 
large-enough power of 2: this happens because one is 
using strongly correlated random numbers to update the 
same lattice sites in successive sweeps (within roughly 
one autocorrelation time) IS S S ■ On the other 
hand, these systematic errors can be eliminated by the 
simple expedient of throwing away one random number 
at the end of each lattice sweep 0, 13, H B S ■ 

It has generally been thought that non-local algorithms 
such as the Swendsen-Wang algorithm ^ and Wolff's 
single-cluster variant ^3 would be immune to these par- 
ticular defects of linear congruential generators, inas- 
much as they employ random numbers in a highly aperi- 
odic way both in "space" and in "time" . We were there- 
fore astonished to find, in our Swendsen-Wang simula- 
tion of the three-dimensional Ising model (13i] , large sys- 
tematic errors on the 128^ and 256'^ lattices that we even- 
tually traced (after much wringing of hands) precisely to 
long-range correlations in the random-number generator. 

Recall that one iteration of the Swendsen-Wang (SW) 
algorithm consists of two steps: first one updates the 
bond occupation variables at a fixed configuration of the 
Ising spin variables; then one computes the connected 
clusters associated to the bond configuration and updates 
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the Ising spin variables by choosing a new spin value in- 
dependently for each cluster. The second (spin-update) 
half of the SW algorithm indeed uses random numbers in 
a thoroughly aperiodic way, because the cluster sizes and 
shapes are random. But the first (bond-update) half uses 
random numbers in a highly structured way: typically 
one sweeps the bonds of the lattice in some simple fixed 
order (e.g. lexicographic). Therefore, if the lattice size 
is very large, the effects of the long-range correlations of 
the random-number generator can be observed within a 
single half-sweep: the random numbers used in updating 
the bonds of one part of the lattice will be strongly corre- 
lated with those used elsewhere in the lattice. One may 
expect this correlation to cause systematic errors partic- 
ularly if (a) the lattice size is commensurate with the lag 
giving rise to long-range correlations (e.g. a power of 2), 
and (b) the system's correlation length is large enough 
so that the long-range correlations of the random-number 
generator couple correlated parts of the lattice. 

The purpose of this note is, first of all, to provide ev- 
idence that such systematic errors can indeed occur and 
that we have accurately diagnosed their origin; and sec- 
ondly, to show how the implementation of the Swendsen- 
Wang algorithm can be modified so as to eliminate (or at 
least radically reduce) these systematic errors. A more 
detailed account will be published elsewhere [T^. 



II. EVIDENCE OF SYSTEMATIC ERRORS 

We simulated the nearest-neighbor three-dimensional 
Ising model on an L x L x L simple-cubic lattice with 
periodic boundary conditions, using the Swendsen-Wang 
(SW) algorithm jlj. We studied lattice sizes L = 4, 6, 8, 
12, 16, 24, 32, 48, 64, 96, 128, 192, 256 and performed be- 
tween 10^ and 10^ SW iterations for each lattice size. We 
did all our runs at /? = 0.22165459, which is Blote et aZ.'s 
best estimate of the critical temperature and is very 
near to the estimates by other workers |1^ J^] (see also 
the review [31 ) ■ We measured a large number of observ- 
ables, including the susceptibility the second-moment 
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L 




deviation (%) 


deviation (cr) 


4 


0.63000(9) 


-0.566% 


-17.48<T 


6 


0.63576(10) 


-0.076% 


3.24<T 


8 


0.63769(10) 


0.004% 


0.44a 


12 


0.63909(11) 


-0.013% 


-0.83cr 


16 


0.63998(8) 


0.002% 


0.15cr 


24 


0.64093(13) 


0.015% 


0.83O- 


32 


0.64122(13) 


-0.010% 


-0.51(7 


48 


0.64172(10) 


-0.007% 


-0.51cr 


64 


0.64223(16) 


0.032% 


1.44cr 


96 


0.64219(16) 


-0.017% 


-0.76O- 


128 


0.62215(25) 


-3.159% 


-78.940- 


192 


0.64383(38) 


0.192% 


3.230- 


256 


0.77798(79) 


21.052% 


169.69(T 


oo 


0.64299(8) 





Table I: Results of our Swendsen-Wang simulations on the 
three-dimensional Ising model at criticality, using a linear con- 
gruential generator with modulus m = 2*®. Error bar (one 
standard deviation) is shown in parentheses. The row marked 
L = oo indicates our best estimate of the asymptotic value 
X* . The last two columns indicate the deviation of each point 
from the fit curve @, in percent and in standard deviations. 
Points deviating by more than 3cr are marked in boldface. 



correlation length ^, the energy E, and the specific heat 
Ch- 

In the first version of our program, the random num- 
bers were supplied by a linear congruential generator 
with modulus m = 2*^, increment c = 1, and multi- 
plier a = 31167285, 10430376854301, 77596615844045 or 
181465474592829. All these multipHers give good results 
on the spectral test in low dimensions, compared to other 
multipliers for the same modulus We verified that 

the runs with the four different multipliers gave results 
that are consistent within error bars for all the major ob- 
servables; after making this verification, we averaged all 
the runs for each L. 

The results for the correlation length ^ are reported 
in the first two columns of Table Finite-size-scaling 
theory predicts that ^/L should behave for large L (if 
indeed we are at the critical temperature) as 

^/L = X* + AL-^ + ... , (1) 

where x* is a universal amplitude ratio characteristic of 
the given system with periodic boundary conditions, lj 
is a correction-to-scaling exponent, A is a nonuniver- 
sal correction-to-scaling amplitude, and the dots indicate 
higher-order corrections to scaling. The data in Tabled 
are qualitatively consistent with (^), except for the points 
at L = 128 and L — 256, which show extremely large de- 
viations. 

A closer examination of the data in Table reveals 
that the point at L = 192 may also exhibit a small but 



statistically significant deviation from the fitting curve. 
To make all these observations more quantitative, let us 
perform a weighted least-squares fit to ^ with uj = 0.82 
(the best estimate from [I3), using all the data with 
Lmin < i < 96 and varying Lmin while checking the 
goodness of fit. A good fit (x^ — 3.85, 6 DF, confidence 
level = 70%) can be obtained already with Lmin = 8, 
yielding 

^/L « 0.64299(8) - 0.02931(79)L"°-^^ (2) 

Not surprisingly, the points L = A and L = 6 show signif- 
icant deviations from the fit curve, due to higher-order 
corrections to scaling. More surprising are the points 
L = 128, which lies roughly 3% (« 79 standard devi- 
ations) below the fit curve, and L — 256, which lies a 
whopping 21% (« 170 standard deviations) above the fit 
curve. Obviously something has gone badly wrong! Fi- 
nally, the point L — 192 lies approximately 0.2% (« 3 
standard deviations) above the fit curve: this may indi- 
cate the presence of a small systematic error also for this 
lattice. 

At first we worried whether we had made a program- 
ming error that might lead to incorrect results on large 
lattices (e.g. due to integer overfiow). We checked the 
program carefully and were unable to find any such mis- 
takes. Moreover, the fact that the systematic discrepancy 
is much smaller (if it exists at all) at L = 192 than at 
L = 128 suggests that the problem — whatever its cause 

— does not arise solely from the lattice being large. 
Intrigued by the fact that these large discrepancies 

might be arising only at lattice sizes that are large pow- 
ers of 2 (or perhaps multiples of large powers of 2), we 
made shorter runs (between 3 x 10'' and 10*^ SW iter- 
ations) at many other lattice sizes — all multiples of 2 
from 4 through 140, and all multiples of 10 through 250 

— in order to check whether any other deviant points 
could be found. The upshot is that — to within the 
statistical error of these shorter runs, which ranges from 
0.2% on small lattices to 1% at L « 128 to an admittedly 
rather crude 2.3% on the largest lattices — there are no 
detectable discrepancies except at L = 128 and 256. 

At L = 128 and L = 256 we found discrepancies not 
only for the correlation length but also for the suscep- 
tibility, the energy and the specific heat. It is a curi- 
ous fact, however, that all the Fortuin-Kasteleyn iden- 
tities [23, equations (3.20)-(3.23)] are verified perfectly 
(to within statistical error). This contrasts with the sys- 
tematic errors found by Damgaard and Heller 8] in a 
Metropolis Monte Carlo simulation of the U{1) Higgs 
model, where a Ward identity was violated by up to 
10 standard deviations, and those found by Ballesteros 
and Martin-Mayor j21J in a Wolff single-cluster simula- 
tion of the two- and three-dimensional Ising models, in 
which Schwinger-Dyson identities were violated by up to 
8 standard deviations. 

For several weeks we had no idea what might be caus- 
ing the systematic discrepancies at L = 128 and L = 256. 
(We felt like a detective with a corpse on his hands but no 
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suspect and no modus operandi.) But then we realized, 
as explained in the Introduction, that long-range corre- 
lations in the random-number generator could cause un- 
desired correlations within a single bond-update sweep. 



III. VARIANT SIMULATIONS 

In order to test whether our proposed explanation for 
the systematic errors is the correct one, we ran vari- 
ant simulations in which two aspects of the simulation 
were systematically altered: the modulus to = 2^ of the 
random-number generator (/3 16, 20, 24, 28, 32, 40, 48, 
60, 63, 64), and the manner in which the random num- 
bers are used within the bond-update subroutine. The 
latter test is essential if we are to prove not only that the 
trouble comes from the random-number generator, but 
more specifically that it comes from the way that the 
random numbers are used in the bond-update subroutine. 

All of the multipliers used here give good results on 
the spectral test in low dimensions compared to other 
multipliers for the same modulus.^ The purpose of try- 
ing random-number generators with less than 48 bits was 
to induce systematic errors on small lattices where they 
could be studied quantitatively to high precision and 
compared with those observed with the 48-bit generator 
on larger lattices. The purpose of trying random-number 
generators with 60/63/64 bits was, of course, to provide 
a standard of comparison in which the systematic error 
is eliminated or at least radically reduced. 

We also tried three variants of the bond-update sub- 
routine: 

Standard: This is our original program, in which the 
bonds are updated in lexicographic order, and one ran- 
dom number is used per bond. 

Aperiodic: Here the bonds are again updated in lexi- 
cographic order, but a random number is used only if the 
two spins are equal. (If the two spins are unequal, the 
corresponding bond is automatically left unoccupied, so 
no random number is needed.) If our explanation of the 
cause of the systematic errors is correct, this strategem 
should eliminate the systematic errors on lattices that 
are multiples of large powers of 2, though it may con- 
ceivably shift those systematic errors to other lattice sizes 
(namely, those for which the lattice size, multiplied by the 
fraction of nearest-neighbor spins that are equal, yields 
a suitable "resonance"). 

Shujfle: The bonds are updated in a random order. ^ 



If our explanation of the cause of the systematic errors 
is correct, this strategem should entirely eliminate the 
systematic errors, even with a relatively poor (e.g. 32- 
bit) random-number generator. 

Our first version of the "shuf&e" subroutine permuted 
the array containing the bond indices. Unfortunately, 
this program ran very slowly — about a factor of 2 slower 
than the "standard" version at L = 16, growing to a 
factor 8 at L — 256 — probably because the highly 
nonlocal access to the bond array caused a large number 
of cache misses. Our second version permuted instead the 
array of random numbers'^; this is statistically equivalent 
but allows the bond array to be accessed in sequential 
order. This program ran less slowly: once again about a 
factor of 2 slower than the "standard" version at L = 16, 
but growing only to a factor « 4 at L = 256. 

The results of all these variant simulations, carried out 
on lattice sizes L = 8^16, 32, 64, 96, 128, 192, 256, will be 
reported elsewhere [lj| ; here we provide only a brief sum- 
mary. We find that the 60/63/64-bit generators give 
consistent results (within statistical error) for all three 
variants of the bond-update subroutine, confirming our 
expectation that they exhibit negligible systematic error 
on lattices L < 256.'* By contrast, each "standard" al- 
gorithm with < 48 bits exhibits detectable systematic 
errors whenever the lattice size L is a multiple of a suffi- 
ciently large power of 2; how large depends on the mod- 
ulus. More precisely, the 16-bit (resp. 20-bit, 24-bit, 28- 
bit, 32-bit, 40-bit, 48-bit) standard algorithm exhibits 
detectable systematic errors whenever L is a multiple of 
8 (resp. 8, 8, 16, 32, 64, 128). In addition, the 48-bit 
"standard" algorithm at L = 192 shows a discrepancy of 
almost 3(J, which may indicate a systematic error. No 
other statistically significant discrepancies are observed. 

We conclude that, if one wants to use a linear congru- 
ential generator with the Swendsen-Wang algorithm, the 
safest approach is to use a generator of 64 bits (or more) 
together with the "shuffle" bond update. Unfortunately, 
the shuffle method is somewhat slow. A much faster — 
and, as far as we can tell, also safe — method is to use 
a 64-bit generator together with the "aperiodic" bond 
update. 

Despite the known problems of linear congruential gen- 
erators arising from long-range correlations, there are still 
several advantages in using them. First, they are rela- 
tively cheap in terms of CPU time, and are convenient 
for use in a series of successive runs because the complete 



1 For moduli 2^2, 2''°, 2"^, 2^°, 2^^ and 2^*, these multipliers 
can be found in published tables of multipliers that perform well 
on the spectral test U jf^j . We double-checked these computa- 
tions, and performed analogous computations from scratch for 
the smaller moduli. 

2 A uniform random permutation of n elements can easily be en- 
erated, in a time of order n, using n — 1 random numbers 
pp. 139-140]. 



^ More precisely, it permuted a LOGICAL array containing the 
results of the comparisons of the random numbers against p = 1 — 
g— 2/3 Tijig requires only one byte storage per bond, rather than 
8 bytes for storing the random number itself, thereby reducing 
both memory usage and cache misses during the generation of 
the random permutation. 

We are continuing runs on the lattice L = 256 in an effort to de- 
tect very small systematic errors. These results will be reported 
later Qj. 
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state of the generator can be saved in a single computer 
word. More importantly, they are well understood theo- 
retically, as re gard s both short-range and long-range 
0, 0) 0j S 1^ ll^l correlations; in particular, excellent 
equidistribution of t-tuples of successive random num- 
bers for small t can be achieved by careful choice of the 
multiplier. By contrast, for more exotic random-number 
generators (e.g. combination generators), the problems 
may not be absent, but simply hidden. 

A more detailed a naly sis of these simulations will be 
published elsewhere along with a discussion of the 
advantages and disadvantages of linear congruential ver- 
sus other types of pseudo-random- number generators (see 

also mill mil). 
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